1 Modelo de Regressão Múltipla


1.1 Estimação dos betas por MQO


O modelo de regressão com \(k>2\) variáveis é conhecido como modelo de regressão múltipla. O modelo de regressão múltipla é uma extensão do modelo de regressão simples. Contudo, as fórmulas matemáticas usadas para calcular \(\hat \beta_1\) e \(\hat \beta_2\) no modelo simples nao podem ser usadas no modelo múltiplo.

Uma forma fácil de se trabalhar com modelos de regressão múltipla é usando algebra matricial. A maior vantagem é que a solução para se estimar os coeficientes pode ser usado para simples e múltipla com qualquer número de variáveis explicativas.

Para se obter a solução do estimador de \(\beta\) na forma matricial, inicialmente se escreve a função de regressão com k-variáveis:

\[ Y_i=\hat \beta_1+ \hat \beta_2 X_{2i}+\hat \beta_3 X_{3i}+ \dots+\hat \beta_k X_{ki}+ \hat u_i \]

que pode ser escrito na forma matricial: com Y sendo um vetor coluna nx1; X uma matriz nxk; \(\hat \beta\) um vetor coluna kx1 e o termo estocástico um vetor coluna nx1.

Os estimadores de MQO são obtidos pela minimização de:

\[ \sum \hat u_i^2= \sum(Y_i- \hat \beta_1 - \hat \beta_2 X_{2i}- \dots -\hat \beta_k X_{ki})^2 \tag{3} \]

com \(\sum \hat u_i^2\) sendo a soma de quadrados dos resíduos. Matricialmente, a soma do quadrados dos resíduos é dada por \(\hat u' \hat u\). A partir de (2), obtem-se

\[ \hat u=Y - X\hat \beta \]

Assim:

\[ \sum \hat u_i^2 = \hat u'\hat u=(Y-X\hat \beta)'(Y-X\hat \beta) \]

\[ =Y'Y-Y'X\hat \beta-\hat \beta' X'Y+ \hat \beta' X'X \hat \beta \]

dado que \(Y'X\hat \beta\) é um escalar, é igual a sua transposta \(\hat \beta'X'Y\). O último termo é a forma quadrática dos elementos de \(\beta\). Então,

\[ \hat u'\hat u=Y'Y-2\hat \beta'X'Y+\hat \beta' X'X \hat \beta \]

que é a função que desejamos minimizar. Para encontrar o ponto de ótimo, deve-se derivar a função e igualá-la a zero.

\[ \frac {\partial \hat u'\hat u}{\partial \hat \beta}=-2X'Y+2X'X\hat \beta=0 \tag{4} \]

\[ X'X\hat \beta=X'Y \tag{5} \]

multiplicando os 2 lados por \((X'X)^{-1}\)

\[ (X'X)^{-1}X'X\hat \beta=(X'X)^{-1}X'Y \tag{6} \]

tem-se

\[ \hat \beta=(X'X)^{-1}X'Y \tag{7} \]


1.2 Demonstração no R - Modelos de Regressão Múltipla

library(wooldridge)
# Carrega a base 'wage2'
data(wage2)
# Remove os valores ausentes (NAs)
sal <- na.omit(wage2)

# Dimensão da base (# linhas  # colunas)
dim(sal)
## [1] 663  17
# Descrição da base
str(sal)
## 'data.frame':    663 obs. of  17 variables:
##  $ wage   : int  769 825 650 562 600 1154 1000 930 1318 1792 ...
##  $ hours  : int  40 40 40 40 40 45 40 43 38 40 ...
##  $ IQ     : int  93 108 96 74 91 111 95 132 119 118 ...
##  $ KWW    : int  35 46 32 27 24 37 44 44 24 47 ...
##  $ educ   : int  12 14 12 11 10 15 12 18 16 16 ...
##  $ exper  : int  11 11 13 14 13 13 16 8 7 9 ...
##  $ tenure : int  2 9 7 5 0 1 16 13 2 9 ...
##  $ age    : int  31 33 32 34 30 36 36 38 28 34 ...
##  $ married: int  1 1 1 1 0 1 1 1 1 1 ...
##  $ black  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ south  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ urban  : int  1 1 1 1 1 0 1 0 1 1 ...
##  $ sibs   : int  1 1 4 10 1 2 1 1 3 1 ...
##  $ brthord: int  2 2 3 6 2 3 1 1 1 1 ...
##  $ meduc  : int  8 14 12 6 8 14 12 13 10 12 ...
##  $ feduc  : int  8 14 12 11 8 5 11 14 10 12 ...
##  $ lwage  : num  6.65 6.72 6.48 6.33 6.4 ...
##  - attr(*, "time.stamp")= chr "25 Jun 2011 23:03"
##  - attr(*, "na.action")= 'omit' Named int [1:272] 2 6 8 12 13 19 20 21 31 36 ...
##   ..- attr(*, "names")= chr [1:272] "2" "6" "8" "12" ...
attach(sal)

O Modelo proposto para ser estimado é:

\[ log(wage)= \beta_0 +\beta_1 educ + \beta_2 exper + \beta_3exper^2+ \beta_4 tenure + u \] onde:

lwage = o logaritmo natural do salário

educ = anos de educação

exper = anos de experiência (trabalhando)

tenure = anos trabalhando com o empregador atual


#Pacote 
library(tseries)

#Estimando a 1 Regressão Multipla Modelo Irrestrito

reg1 <- lm(log(wage) ~ educ + exper + I(exper^2) + tenure, data=sal)
summary(reg1)
## 
## Call:
## lm(formula = log(wage) ~ educ + exper + I(exper^2) + tenure, 
##     data = sal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8515 -0.2351  0.0170  0.2491  1.3157 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.5430834  0.1419558  39.048  < 2e-16 ***
## educ        0.0763886  0.0075169  10.162  < 2e-16 ***
## exper       0.0079634  0.0154608   0.515  0.60667    
## I(exper^2)  0.0004488  0.0006674   0.673  0.50149    
## tenure      0.0095632  0.0030992   3.086  0.00212 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3791 on 658 degrees of freedom
## Multiple R-squared:  0.1591, Adjusted R-squared:  0.154 
## F-statistic: 31.13 on 4 and 658 DF,  p-value: < 2.2e-16
#Testes de Normalidade
jarque.bera.test(reg1$residuals)
## 
##  Jarque Bera Test
## 
## data:  reg1$residuals
## X-squared = 25.853, df = 2, p-value = 2.433e-06
#Estimando a 1 Regressão Multipla Modelo Restrito

regrest1 <- lm(log(wage) ~ educ  + tenure, data=sal)
summary(regrest1)
## 
## Call:
## lm(formula = log(wage) ~ educ + tenure, data = sal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.90839 -0.24389  0.02665  0.25641  1.30413 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.884887   0.095844  61.401  < 2e-16 ***
## educ        0.060522   0.006696   9.039  < 2e-16 ***
## tenure      0.014057   0.002955   4.757 2.42e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3843 on 660 degrees of freedom
## Multiple R-squared:  0.1337, Adjusted R-squared:  0.131 
## F-statistic: 50.91 on 2 and 660 DF,  p-value: < 2.2e-16
#Regressão com variáveis em log
reg2<-lm(log(wage) ~ educ + log(exper) + tenure, data=sal)
summary(reg2)
## 
## Call:
## lm(formula = log(wage) ~ educ + log(exper) + tenure, data = sal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.88595 -0.23687  0.02024  0.24902  1.29614 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.500423   0.145736  37.742  < 2e-16 ***
## educ        0.069977   0.007174   9.754  < 2e-16 ***
## log(exper)  0.116374   0.033440   3.480 0.000534 ***
## tenure      0.011679   0.003009   3.881 0.000115 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3811 on 659 degrees of freedom
## Multiple R-squared:  0.1493, Adjusted R-squared:  0.1454 
## F-statistic: 38.55 on 3 and 659 DF,  p-value: < 2.2e-16
jarque.bera.test(reg2$residuals)
## 
##  Jarque Bera Test
## 
## data:  reg2$residuals
## X-squared = 33.743, df = 2, p-value = 4.707e-08
data("econmath")
summary(econmath)
##       age             work            study           econhs      
##  Min.   :18.00   Min.   : 0.000   Min.   : 0.00   Min.   :0.0000  
##  1st Qu.:19.00   1st Qu.: 0.000   1st Qu.: 8.50   1st Qu.:0.0000  
##  Median :19.00   Median : 8.000   Median :12.00   Median :0.0000  
##  Mean   :19.41   Mean   : 8.626   Mean   :13.92   Mean   :0.3703  
##  3rd Qu.:20.00   3rd Qu.:15.000   3rd Qu.:18.00   3rd Qu.:1.0000  
##  Max.   :29.00   Max.   :37.500   Max.   :50.00   Max.   :1.0000  
##                                                                   
##      colgpa          hsgpa           acteng          actmth     
##  Min.   :0.875   Min.   :2.355   Min.   :12.00   Min.   :12.00  
##  1st Qu.:2.446   1st Qu.:3.110   1st Qu.:20.00   1st Qu.:20.00  
##  Median :2.813   Median :3.321   Median :23.00   Median :23.00  
##  Mean   :2.815   Mean   :3.342   Mean   :22.59   Mean   :23.21  
##  3rd Qu.:3.207   3rd Qu.:3.589   3rd Qu.:25.00   3rd Qu.:26.00  
##  Max.   :4.000   Max.   :4.260   Max.   :34.00   Max.   :36.00  
##                                  NA's   :42      NA's   :42     
##       act           mathscr            male        calculus     
##  Min.   :13.00   Min.   : 0.000   Min.   :0.0   Min.   :0.0000  
##  1st Qu.:21.00   1st Qu.: 7.000   1st Qu.:0.0   1st Qu.:0.0000  
##  Median :23.00   Median : 8.000   Median :0.5   Median :1.0000  
##  Mean   :23.12   Mean   : 7.875   Mean   :0.5   Mean   :0.6764  
##  3rd Qu.:25.00   3rd Qu.: 9.000   3rd Qu.:1.0   3rd Qu.:1.0000  
##  Max.   :33.00   Max.   :10.000   Max.   :1.0   Max.   :1.0000  
##  NA's   :42                                                     
##      attexc          attgood          fathcoll         mothcoll     
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :0.0000   Median :1.0000   Median :1.0000   Median :1.0000  
##  Mean   :0.2967   Mean   :0.5864   Mean   :0.5245   Mean   :0.6285  
##  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##                                                                     
##      score      
##  Min.   :19.53  
##  1st Qu.:64.06  
##  Median :74.22  
##  Mean   :72.60  
##  3rd Qu.:82.79  
##  Max.   :98.44  
## 
#Estimação do modelo 
reg3<-lm(score ~ study + age + actmth, data=econmath)
summary(reg3)
## 
## Call:
## lm(formula = score ~ study + age + actmth, data = econmath)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -45.032  -7.602   1.660   8.904  27.893 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 56.38368    9.21953   6.116 1.49e-09 ***
## study        0.07111    0.05472   1.299   0.1942    
## age         -0.94664    0.45149  -2.097   0.0363 *  
## actmth       1.44843    0.11330  12.784  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.13 on 810 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.1725, Adjusted R-squared:  0.1694 
## F-statistic: 56.27 on 3 and 810 DF,  p-value: < 2.2e-16
#Teste de Normalidade
shapiro.test(reg3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg3$residuals
## W = 0.97894, p-value = 1.919e-09
#Estimação do modelo com variáveis transformadas
reg4<-lm(score~study+log(age)+log(actmth), data=econmath)
summary(reg4)
## 
## Call:
## lm(formula = score ~ study + log(age) + log(actmth), data = econmath)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.877  -7.698   1.570   8.796  27.725 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20.70591   28.87977   0.717   0.4736    
## study         0.07665    0.05462   1.403   0.1609    
## log(age)    -18.23167    9.28094  -1.964   0.0498 *  
## log(actmth)  33.50102    2.58237  12.973   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.1 on 810 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.1765, Adjusted R-squared:  0.1734 
## F-statistic: 57.87 on 3 and 810 DF,  p-value: < 2.2e-16
#Teste de Normalidade
shapiro.test(reg4$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg4$residuals
## W = 0.97933, p-value = 2.536e-09
data("crime2")

#Estimação do modelo 
reg5 <- lm(crimes ~ pop + unem + officers + pcinc, data=crime2)
summary(reg5)
## 
## Call:
## lm(formula = crimes ~ pop + unem + officers + pcinc, data = crime2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -27745  -6110  -1200   4139  60939 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.523e+03  1.058e+04  -0.333    0.740    
## pop          6.432e-02  1.079e-02   5.961 5.23e-08 ***
## unem        -6.381e+01  5.519e+02  -0.116    0.908    
## officers     1.516e+01  3.596e+00   4.216 6.06e-05 ***
## pcinc        4.779e-01  8.211e-01   0.582    0.562    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12600 on 87 degrees of freedom
## Multiple R-squared:  0.828,  Adjusted R-squared:  0.8201 
## F-statistic: 104.7 on 4 and 87 DF,  p-value: < 2.2e-16
#Verificação dos Resíduos
qqnorm(reg5$residuals)
qqline(reg5$residuals)

#Teste de Normalidade
shapiro.test(reg4$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg4$residuals
## W = 0.97933, p-value = 2.536e-09
#Estimação do modelo com variáveis transformadas
reg6<-lm(log(crimes)~log(pop)+log(unem)+log(officers)+log(pcinc), data=crime2)
summary(reg6)
## 
## Call:
## lm(formula = log(crimes) ~ log(pop) + log(unem) + log(officers) + 
##     log(pcinc), data = crime2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.46042 -0.16369 -0.00741  0.16572  0.62182 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.25623    1.61835   0.776    0.440    
## log(pop)       0.52707    0.11930   4.418 2.86e-05 ***
## log(unem)     -0.08666    0.10124  -0.856    0.394    
## log(officers)  0.43781    0.10527   4.159 7.47e-05 ***
## log(pcinc)    -0.03178    0.16067  -0.198    0.844    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2628 on 87 degrees of freedom
## Multiple R-squared:  0.8801, Adjusted R-squared:  0.8745 
## F-statistic: 159.6 on 4 and 87 DF,  p-value: < 2.2e-16
#Verificação dos Resíduos
qqnorm(reg6$residuals)
qqline(reg6$residuals)

#Teste de Normalidade
shapiro.test(reg6$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg6$residuals
## W = 0.98173, p-value = 0.2253

2 Quebra de Pressupostos do Modelo Clássico de Regressão Linear


2.1 Multicolinearidade


Uma das premissas do Modelo Clássico de Regressão Linear é que não existe multicolinearidade entre as variáveis explicativas no modelo de regressão a ser estimado. Originalmente, multicolinearidade designava a existência de uma “relação perfeita” entre algumas ou todas as variáveis explicativas de um modelo de regressão. Atualmente, também é considerado o caso de multicolinearidade menos que perfeita.


x1 <- c(10,15,18,24,30)
x2 <- c(50,75,90,120,150)
x3 <- c(52,75,97,129,152)
x <- cbind(x1,x2,x3)
cor(x)
##           x1        x2        x3
## x1 1.0000000 1.0000000 0.9958817
## x2 1.0000000 1.0000000 0.9958817
## x3 0.9958817 0.9958817 1.0000000


A correlação entre \(X_1\) e \(X_2\) é igual a 1. A correlação entre \(X_1\) e \(X_3\) é 0,995. No primeiro caso se tem perfeita correlação. Se fosse estimar um modelo de regressão com estas duas variáveis explicativas, \(X_1\) e \(X_2\), teriamos um caso de multicolinearidade perfeita e os betas não seriam estimáveis. A demonstração é dada em Gujarati & Porter (2009).


A título de exemplo, se pode observar com a matrix X abaixo:


\[ \mathbf{X}=\left[\begin{array}{ccc} 1 & 10 & 50 \\ 1 & 15 & 75 \\ 1 & 18 & 90 \\ \end{array} \right]; \mathbf{X'X}=\left[\begin{array}{ccc} 3 & 43 & 215 \\ 43 & 649 & 3245 \\ 215 & 3245 & 16225 \\ \end{array} \right] \]


sendo o determinante desta matriz igual a zero \(|\mathbf{X'X}|=0\).


x1 <- c(1,1,1)
x2 <- c(10,15,18)
x3 <- c(50,75,90)
x <- cbind(x1,x2,x3)

x <- as.matrix(x)

# X'X
m1 <- t(x) %*% x

m1
##     x1   x2    x3
## x1   3   43   215
## x2  43  649  3245
## x3 215 3245 16225
# Determinante

det(m1)
## [1] 0


A colinearidade se refere às relações lineares entre as variáveis explicativas. Não inclui relações não lineares, como por exemplo:


\[ Y_i=\beta_0+\beta_1X+\beta_2X^2+\beta_3X^3 \]


em que Y é o custo total e X a quantidade produzida


2.2 Consequências da multicolinearidade


Os estimadores de MQO continuam não viesados e eficientes (possuem variância mínima).

Contudo, esta variância é “inflada”, ou seja, muito grande.

Se o erro padrão é grande, os intervalos de confiança também serão mais amplos:


\[ IC=\hat\beta\pm t_{\alpha/2}ep(\hat \beta) \]


Como o erro padrão é grande, normalmente os betas são não significativos, ou seja, iguais a zero.

O \(R^2\) tende a ser muito alto.

Os estimadores de MQO e os erros-padrões podem ser sensíveis às pequenas mudanças nos dados.

Em resumo, os resultados encontrados se tornam duvidosos quando se regride com regressores colineares.


2.3 Como detectar a multicolinearidade?


Na regressão, encontrar muitos betas estimados não significativos e o \(R^2\) muito alto.

Altos valores de correlações entre pares de regressores;

Calcular regressões auxiliares. Dado o modelo:


\[ Y_i=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_3+u \]


calcular o \(R^2\) e chamar de \(R^2\) geral. Depois regredir:


\[ X_1=\beta_0+\beta_2X_2+\beta_3X_3+u \]

\[ X_2=\beta_0+\beta_1X_1+\beta_3X_3+u \]

\[ X_3=\beta_0+\beta_1X_1+\beta_2X_2+u \]

calcular em cada regressão auxiliar os \(R^2\): a) \(R^2_{x_2.x_3x_4}\); b) \(R^2_{x_3.x_2x_4}\); c) \(R^2_{x_4.x_2x_3}\)


comparar cada um com o \(R^2\) geral. Pela Regra Prática de Klein, a multicolinearidade é um problema sério se o \(R^2\) obtido em todas as regressões auxiliares for maior que o \(R^2\) geral.


Fator de Inflação de Variância - Considerando uma regressão múltipla com duas variáveis explicativas, as variâncias dos coeficientes de inclinação são dadas por:


\[ Var(\hat\beta_2)=\frac{\sigma^2}{\sum x^2_2(1-r^2_{23})} \]

\[ Var(\hat\beta_3)=\frac{\sigma^2}{\sum x^2_3(1-r^2_{23})} \]


em que \(\sigma^2\) é a variância do termo de erro e \(r^2_{23}\) é o coeficiente de correlação entre \(x_2\) e \(x_3\).

se considerarmos a razão \(\frac{1}{1-r^2_{23}}\), esta mensura o grau na qual a variância do estimador de MQO é inflado devido a colinearidade. Este fator é conhecido como FIV (Fator de Inflação de Variância).

Se o FIV for maior do que 10 diz-se que esta variável é altamente colinear; Com mais de 2 variáveis, o quadrado do coeficiente de correlação pode ser obtido do \(R^2\) das regressões auxiliares.


2.4 Medidas corretivas para Multicolinearidade


  1. Aumentar o tamanho da amostra;

  2. Exclusão de variáveis, desde que não cause viés de especificação ou omissão de variável relevante.


2.5 Multicolinearidade: Exemplo no R


2.5.1 Estimação do Modelo

#Inicio do Script
#Pacotes a serem utilizados
library(car)
library(corrplot)
library(HH)
library(lmtest)
library(sandwich)

#Regressao Multipla
#tenure= anos no emprego atual
regressao1 <- lm(log(wage) ~ educ + log(exper) + tenure, data=sal)
summary(regressao1)
## 
## Call:
## lm(formula = log(wage) ~ educ + log(exper) + tenure, data = sal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.88595 -0.23687  0.02024  0.24902  1.29614 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.500423   0.145736  37.742  < 2e-16 ***
## educ        0.069977   0.007174   9.754  < 2e-16 ***
## log(exper)  0.116374   0.033440   3.480 0.000534 ***
## tenure      0.011679   0.003009   3.881 0.000115 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3811 on 659 degrees of freedom
## Multiple R-squared:  0.1493, Adjusted R-squared:  0.1454 
## F-statistic: 38.55 on 3 and 659 DF,  p-value: < 2.2e-16


2.5.2 Início da análise: verificação da correlação

#Matriz de Correlacoes das variaveis explicativas
x <- wage1[, c(2,3,4)]
cor(x)
##               educ      exper      tenure
## educ    1.00000000 -0.2995418 -0.05617257
## exper  -0.29954184  1.0000000  0.49929145
## tenure -0.05617257  0.4992914  1.00000000
corrplot(cor(x), order="hclust", tl.col="black", tl.cex = .75)


2.5.3 Estimação das regressões auxiliares

#Regressoes Auxiliares

aux1 <- lm(educ ~ exper + tenure)
summary(aux1)
## 
## Call:
## lm(formula = educ ~ exper + tenure)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0084 -1.4863 -0.3875  1.5734  5.1828 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.21282    0.22610  71.708   <2e-16 ***
## exper       -0.25293    0.01889 -13.389   <2e-16 ***
## tenure       0.04850    0.01591   3.048   0.0024 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.981 on 660 degrees of freedom
## Multiple R-squared:  0.2143, Adjusted R-squared:  0.2119 
## F-statistic: 90.01 on 2 and 660 DF,  p-value: < 2.2e-16
aux2 <- lm(exper~ educ + tenure)
summary(aux2)
## 
## Call:
## lm(formula = exper ~ educ + tenure)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.2109  -2.2600   0.0697   2.5506   9.8563 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 21.26670    0.90280  23.556  < 2e-16 ***
## educ        -0.84445    0.06307 -13.389  < 2e-16 ***
## tenure       0.23308    0.02784   8.373 3.37e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.619 on 660 degrees of freedom
## Multiple R-squared:  0.2797, Adjusted R-squared:  0.2776 
## F-statistic: 128.2 on 2 and 660 DF,  p-value: < 2.2e-16
aux3 <- lm(tenure ~ educ + exper)
summary(aux3)
## 
## Call:
## lm(formula = tenure ~ educ + exper)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.693 -3.985 -0.161  3.671 13.191 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.3929     1.6275  -0.856   0.3924    
## educ          0.2862     0.0939   3.048   0.0024 ** 
## exper         0.4120     0.0492   8.373 3.37e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.812 on 660 degrees of freedom
## Multiple R-squared:  0.09684,    Adjusted R-squared:  0.0941 
## F-statistic: 35.38 on 2 and 660 DF,  p-value: 2.524e-15


2.5.4 Cálculo do Fator de Inflação da Variância

#Fator de Inflacao da Variancia (FIV)
vif(regressao1)
##       educ log(exper)     tenure 
##   1.168469   1.227742   1.055337

2.5.5 Uso do pacote {performance}

pacman::p_load(
"performance",
"lme4",
"see",
"qqplotr"
)
## 
## The downloaded binary packages are in
##  /var/folders/pw/v66yct9s70s49dj4xyx5rk4m0000gn/T//RtmpOYPLUS/downloaded_packages
model_performance(regressao1)
## # Indices of model performance
## 
## AIC     |    R2 | R2 (adj.) |  RMSE | Sigma
## -------------------------------------------
## 608.169 | 0.149 |     0.145 | 0.380 | 0.381
check_model(regressao1)


2.6 Heterocedasticidade


Ocorre quando a variância do termo de erro não é constante \((\sigma^2_i)\). As fontes principais são: presença de outliers, modelo incorretamente especificado, assimetria na distribuição de um ou mais regressores e erros de digitação dos dados.

É mais comum em dados de seção cruzada.

Com as relação às propriedades dos estimadores de MQO, na presença de heterocedasticidade continuam não-viesados e consistentes, mas não são mais eficientes. O estimador eficiente, ou seja, com variância mínima é o de Mínimos Quadrados Generalizados (MQG).


2.7 Consequências da Heterocedasticidade


O uso do método de MQO na presença de heterocedasticidade leva a resultados incorretos das significâncias (Testes t e F), levando a conclusões e inferências que podem estar equivocadas.


2.8 Como detectar a Heterocedasticidade?


Sempre uma boa análise gráfica dos resíduos é importante. Contudo, existem testes formais para análise de heterocedasticidade. Os principais são os de White e Breusch-Pagan-Godfrey (BPG).


2.9 Exemplo no R


#Regressao Multipla
regressao1 <- lm(log(wage) ~ urban + married
                 + educ + exper, data=sal)
summary(regressao1)
## 
## Call:
## lm(formula = log(wage) ~ urban + married + educ + exper, data = sal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.97210 -0.20965  0.02298  0.22976  1.19926 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.217673   0.131350  39.723  < 2e-16 ***
## urban       0.198629   0.031859   6.235 8.09e-10 ***
## married     0.218702   0.047801   4.575 5.69e-06 ***
## educ        0.074823   0.007182  10.418  < 2e-16 ***
## exper       0.020461   0.003760   5.442 7.44e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3665 on 658 degrees of freedom
## Multiple R-squared:  0.2144, Adjusted R-squared:  0.2097 
## F-statistic: 44.91 on 4 and 658 DF,  p-value: < 2.2e-16
#Teste de Normalidade dos residuos
jarque.bera.test(regressao1$residuals)
## 
##  Jarque Bera Test
## 
## data:  regressao1$residuals
## X-squared = 72.811, df = 2, p-value < 2.2e-16
shapiro.test(regressao1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  regressao1$residuals
## W = 0.98616, p-value = 6.381e-06
#Analise de Multicolinearidade
#Matriz de Correlacoes das variaveis explicativas
x <- sal[, c(5,6,9,12)]
cor(x)
##                educ       exper     married       urban
## educ     1.00000000 -0.45082581 -0.05671975  0.10016540
## exper   -0.45082581  1.00000000  0.09847967 -0.03804240
## married -0.05671975  0.09847967  1.00000000 -0.03942462
## urban    0.10016540 -0.03804240 -0.03942462  1.00000000
corrplot(cor(x), order="hclust", tl.col="black", tl.cex = .75)

#Calculo do Fator de Inflacao da Variancia
vif(lm(log(wage) ~ urban + married + educ + exper, data=sal))
##    urban  married     educ    exper 
## 1.011420 1.011208 1.266198 1.263695


2.9.1 Verificação dos resíduos


#Analise dos residuos do modelo
par(mfrow=c(2,2))
plot(regressao1)

#Salvando os residuos do modelo e os valores estimados
resid_sq <- (regressao1$residuals)^2
ajustados1 <- regressao1$fitted.values
par(mfrow=c(1,1))
plot(ajustados1, resid_sq)


2.9.2 Heterocedasticidade: teste de White


A hipótese nula do Teste de White é os resíduos são homocedásticos. Após a regressão, pegar os resíduos ao quadrado e estimar contra os regressores e estes ao quadrado. é possível incluir produtos cruzados também. A estatística de teste \(nR^2\) segue uma distribuição de qui-quadrado com graus de liberdade igual ao numero de parâmetros estimados. É um teste para grandes amostras.


#Testes de Heterocedasticidade de White
#Caso especial do BPG
#H0: Residuos sao homocedasticos.
regres_white <- lm(log(wage) ~ urban + married
                 + educ + exper + I(educ^2) + I(exper^2), data=sal)
bptest(regres_white) 
## 
##  studentized Breusch-Pagan test
## 
## data:  regres_white
## BP = 10.247, df = 6, p-value = 0.1146


2.9.3 Heterocedasticidade: teste de BPG


Semelhante ao Teste de White, mas os resíduos ao quadrado são estimados no teste BPG contra os regressores originais apenas. Na estatística de teste é possível usar o valor do F ou então o \(R^2\) multiplicado pelo tamanho da amostra “n”. A estatística de teste segue a distribuição de qui-quadrado com graus de liberdade igual ao número de regressores no modelo. A hipótese nula é homocedasticidade.


#Testes de Heterocedasticidade de White
#Caso especial do BPG
#H0: Residuos sao homocedasticos.
#Testes de Heterocedasticidade de Breusch-Pagan
bptest(regressao1)
## 
##  studentized Breusch-Pagan test
## 
## data:  regressao1
## BP = 9.7202, df = 4, p-value = 0.04541
#Testes de Heterocedasticidade NCV
ncvTest(regressao1)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 3.358497, Df = 1, p = 0.066859


2.9.4 Uso do pacote {performance}


pacman::p_load(
"performance",
"lme4",
"see",
"qqplotr"
)
## 
## The downloaded binary packages are in
##  /var/folders/pw/v66yct9s70s49dj4xyx5rk4m0000gn/T//RtmpOYPLUS/downloaded_packages
model_performance(regressao1)
## # Indices of model performance
## 
## AIC     |    R2 | R2 (adj.) |  RMSE | Sigma
## -------------------------------------------
## 557.345 | 0.214 |     0.210 | 0.365 | 0.366
check_model(regressao1)


2.10 Heterocedasticidade: como corrigir?


Usar mínimos quadrados ponderados se conhecer a verdadeira variância heterocedástica \(\sigma^2_i\). Contudo, raramente isto ocorre. O que mais se usa é transformação logarítmica de variáveis com grande variância e o procedimento de White para estimar erros-padrão robustos. O procedimento não altera a estimação no ponto, apenas os erros-padrão corrigindo a heterocedasticidade.

A correção clássica de White para a matriz de variância dos coeficientes é dada por:

\[ Var(\hat \beta|X)=(X'X)^{-1}X'diag(e^2)X(X'X)^{-1} \] em que \(e^2\) são os resíduos ao quadrado e \(X\) é a matriz de variáveis explicativas do modelo.

Contudo, existem diversos outros métodos de ajustamento desta fórmula para a matriz de var-cov dos resíduos, denotada por \(\Omega\). No R, a opção “const” retorna os resultados de \(\sigma^2(X'X)^{-1}\). As outras opções retornam estimadores de White na forma da expressão e diferem em termos de \((X'X)^{-1}X' \Omega X(X'X)^{-1}\) e diferem em como é \(\Omega\). A diagonal de \(\Omega\) será formada pelos elementos \(\omega_i\) exposto para cada opção de HC (hc0, hc1, hc2, hc3, hc4):


  1. const = \(\omega_i= \hat \sigma^2\)
  2. HC0 = \(\omega_i= \hat u_i^2\) é matriz clássica de correção de Eicker (1963) e popularizada por White (1980);
  3. HC1 = \(\omega_i= \frac {n}{n-k}\hat u_i^2\)
  4. HC2 = \(\omega_i= \frac {\hat u_i^2}{1-h_i}\)
  5. HC3 = \(\omega_i= \frac {\hat u_i^2}{(1-h_i)^2}\)
  6. HC4 = \(\omega_i= \frac {\hat u_i^2}{(1-h_i)^{\delta_i}}\) é a matriz de correção conforme Cribari-Neto (2004) para aperfeiçoar a performance em pequenas amostras com presença de observações influentes.

“hc1”, “hc2” e “hc3” as matrizes de correção sugeridas por MacKinnon e White (1985) e aperfeiçoadas para pequenas amostras conforme Long e Ervin (2000).

para \(h_i=H_{ii}\) como os elementos diagonais da matriz estimada, \(\bar h\) é sua média é \(\delta_i=min {\{4, h_i/ \bar h}\}\). A opção que retorna os mesmos resultados que o “padrão” de correção de White no Stata e Eviews é a “hc1”.


#Correcao da Heterocedasticidade
regressao1 <- lm(log(wage) ~ urban + married
                 + educ + exper, data=sal)
coeftest(regressao1, vcov=vcovHC, type='HC1') #(Eviews)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 5.2176734  0.1338840 38.9716 < 2.2e-16 ***
## urban       0.1986294  0.0313196  6.3420 4.216e-10 ***
## married     0.2187024  0.0507873  4.3062 1.914e-05 ***
## educ        0.0748232  0.0072903 10.2635 < 2.2e-16 ***
## exper       0.0204610  0.0037196  5.5008 5.417e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(regressao1, vcov=vcovHC, type='HC4') #sugerido (2011)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 5.2176734  0.1346921 38.7378 < 2.2e-16 ***
## urban       0.1986294  0.0314178  6.3222 4.759e-10 ***
## married     0.2187024  0.0517175  4.2288 2.682e-05 ***
## educ        0.0748232  0.0073184 10.2240 < 2.2e-16 ***
## exper       0.0204610  0.0037343  5.4792 6.089e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


3 Autocorrelação e viés de especificação


3.1 O que é Autocorrelação?

Pode existir autocorrelação no cross-section, mas é muito mais presente em séries temporais. Existe dependência entre os termos de erro \((E(u_i u_j \ne 0 \quad \forall \quad i \ne j))\). Diversos fatores podem gerar a correlação serial, como por exemplo:

  1. inércia;

  2. viés de especificação;

  3. ausência de estacionariedade;

Com relação à propriedade dos estimadores, com autocorrelação, os estimadores de MQO são não-viesados e consistentes, mas não são eficientes. Os testes t e F não são mais válidos.

Nesta situação, da mesma forma que para a heterocedasticidade, é melhor utilizar o método de Mínimos Quadrados Generalizados (MQG).

Considere um modelo com resíduos que tem correlação de primeira ordem:

\[ u_t=\rho u_{t-1}+v_t \] em que \(\rho\) é o parâmetro de autocorrelação e “v” é um termo de erro “bem comportado”, ou seja, não autocorrelacionado que segue a distribuição Normal com média zero e variância constante \(\sigma^2\).

O coeficiente de autocorrelação \(\rho\) pode ser obtido através de

\[ \hat \rho = \frac {cov(\epsilon_t, \epsilon_{t-1})}{[var(\epsilon_t)^{1/2}] [var(\epsilon_{t-1})^{1/2}]} \]

Com base nesta idéia, alguns testes foram desenvolvidos visando verificar a presença de autocorrelação nos resíduos da regressão.

3.2 Autocorrelação: Como Detectar?

Considere o modelo abaixo, uma série temporal que se inicia em 1947 e faz uma regressão do consumo em função da renda, da riqueza e da taxa de juros. As três primeiras variáveis estão em log.

#Inicio do Script
#Pacotes a serem usados
library(forecast)
library(ggplot2)
library(scatterplot3d)

options(digits=4)

data("consump")
dados <- consump
dados <- ts(dados, start=c(1959))

#Estimacao da Equacao de Regressao
regressao1 <- lm(log(c) ~log(y) + inf, data=dados)
summary(regressao1)
## 
## Call:
## lm(formula = log(c) ~ log(y) + inf, data = dados)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.03009 -0.00695  0.00240  0.00847  0.01672 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.264900   0.077930    3.40   0.0017 ** 
## log(y)       0.951550   0.008284  114.87   <2e-16 ***
## inf         -0.001872   0.000648   -2.89   0.0067 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0115 on 34 degrees of freedom
## Multiple R-squared:  0.998,  Adjusted R-squared:  0.998 
## F-statistic: 7.24e+03 on 2 and 34 DF,  p-value: <2e-16


3.2.1 MÉTODO GRÁFICO


Apesar de existirem testes estatísticos, uma análise gráfica pode ser útil para tentar verificar algum comportamento nos dados que de indicativo se existência de autocorrelação.

#Analise dos residuos do modelo
resid <- regressao1$residuals
resid <- ts(resid, start=c(1947))
par(mfrow=c(1,1))
scatterplot3d(lag(resid,-1), resid, pch=16, box=FALSE)


3.2.2 TESTE DE DURBIN-WATSON


O teste mais conhecido é o “d” de Durbin-Watson (razão da soma das diferenças ao quadrado entre os sucessivos resíduos e a soma de quadrados dos resíduos):

\[ d=\frac{\sum_{t=2}^{t=n}(\hat{u}_t-\hat{u}_{t-1})^2}{\sum_{t=1}^{t=n}(\hat{u}_t^2)} \]

Possui importantes premissas:

  1. a regressão precisa ter intercepto;

  2. o termo de erro precisa ser um AR(1) (\(u_t=\rho u_{t-1}+v_t\)), não pode ser AR de ordem superior;

  3. o erro deve ser normal;

  4. o modelo não pode ter termos defasados da variável dependente (y);

Não há um único valor crítico que leve a rejeição ou não da hipótese nula (ausência de autocorrelação), contudo o R mostra um valor de Probabilidade. Gujarati e Porter (2009) demonstram facilmente que \(d\approx2(1- \hat{\rho})\).

\[ \hat{u}_t=\hat{\rho} u_{t-1}+v_t \]

Dado que \(-1\le \rho \le 1\), tem-se que \(0\le d \le 4\), que são os limites de “d”. Como regra prática, espera-se que d esteja próximo de 2 para pressupor que não existe autocorrelação nos resíduos. Quanto mais próximo de 0 (zero), maior a evidência de auto positiva e de 4, de auto negativa.

#Teste de Durbin-Watson
# Testa the hipotese que nao existe correlacao de 1 lag nos residuos
dwtest(regressao1, alt="two.sided")
## 
##  Durbin-Watson test
## 
## data:  regressao1
## DW = 0.82, p-value = 1e-05
## alternative hypothesis: true autocorrelation is not 0


3.2.3 TESTE DE BREUSCH-GODFREY


Teste que permite regressandos defasados inclusos como regressores, processos autoregressivos (AR) de ordens maiores do que 1 e Termos de médias móveis (MA).

Suponha que o termo de erro siga a seguinte estrutura:

\[ u_t=\rho_1u_{t-1}+\rho_2u_{t-2}+\dots+\rho_pu_{t-p}+v_t \]

\(v_t\) é um erro que segue as hipóteses do MCRL.

A hipótese nula do teste BG é

\[ H_0=\rho_1=\rho_2=\dots=\rho_p=0 \]

Ausência de autocorrelação de qualquer ordem. Dada a não observação de \(u_t\) é usado a sua estimativa e estimada uma regressão auxiliar de \(e_t\) contra as variáveis explicativas e os termos de erros defasados. Em grandes amostras, o teste BG segue a distribuição \((n-p)R^2 \sim \chi_p^2\).

O teste BG possui as seguintes características:

  1. A variância de \(u_t\) seja homocedástica;

  2. Definir a ordem de “p”. Normalmente isto é um processo de tentativa e erro, com o uso de critérios de Akaike ou Schwarz para definição do melhor modelo. Deve ser escolhido o modelo com menor valor de critério de informação;


#Teste de Breusch-Godfrey
#AIC(regressao1)
#BIC(regressao1)
bg1 <- bgtest(regressao1,1)
coeftest(bg1)
## 
## z test of coefficients:
## 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.007700   0.063945    0.12     0.90    
## log(y)       -0.000973   0.006798   -0.14     0.89    
## inf           0.000351   0.000538    0.65     0.51    
## lag(resid)_1  0.597970   0.142777    4.19  2.8e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bg1
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  regressao1
## LM test = 13, df = 1, p-value = 3e-04
bg2 <- bgtest(regressao1,2)
coeftest(bg2)
## 
## z test of coefficients:
## 
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.006783   0.064787    0.10   0.9166   
## log(y)       -0.000892   0.006887   -0.13   0.8970   
## inf           0.000391   0.000553    0.71   0.4790   
## lag(resid)_1  0.555076   0.175969    3.15   0.0016 **
## lag(resid)_2  0.078814   0.184311    0.43   0.6689   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bg2
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  regressao1
## LM test = 13, df = 2, p-value = 0.002
bg3 <- bgtest(regressao1,3)
coeftest(bg3)
## 
## z test of coefficients:
## 
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.009099   0.065366    0.14   0.8893   
## log(y)       -0.001105   0.006946   -0.16   0.8736   
## inf           0.000309   0.000569    0.54   0.5863   
## lag(resid)_1  0.563184   0.177685    3.17   0.0015 **
## lag(resid)_2  0.141504   0.205324    0.69   0.4907   
## lag(resid)_3 -0.135297   0.188919   -0.72   0.4739   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bg3
## 
##  Breusch-Godfrey test for serial correlation of order up to 3
## 
## data:  regressao1
## LM test = 13, df = 3, p-value = 0.004


3.3 Autocorrelação: Como Corrigir?


Existem diversas formas de corrigir a autocorrelação:

  1. Transformar as séries, diferenciando-as, caso se saiba o valor de \(\rho\) ou considerando \(\rho=1\);

  2. Fazer uma transformação generalizada: \(Y_t- \hat{\rho}(Y_{t-1})\) para obter \(\hat{\rho}\), fazer \(\approx 1-\frac{d}{2}\);

  3. Usar o método de correção dos erros padrões de Newey-West (HAC-Consistente com auto e hetero), desde que a amostra seja grande. Erros robustos.

Com base em Figueiro (2021), o cálculo dos erros-padrão por estimativas robustas serão desejáveis para formas mais gerais de correlação serial, diga-se, com autocorrelações de ordens superiores.

Desta forma, estima-se o modelo padrão de regressão linear por MQO. Em seguida, se pega os resíduos estimados de uma regressão auxiliar de \(x_{1t}\) em função dos demais \(x_{kt}\) e calcula-se \(\hat a = \hat r_t \hat u_t\) em que \(\hat u_t\) são os resíduos da estimação original por MQO. Por uma escolha de “g”, que pode variar entre a parte inteira de \(4(n/100)^{2/9}\) ou \(n^{1/4}\), calcula-se


\[ \hat v= \sum_{t=1}^{n} k^2 \hat a^2_t + 2 \sum_{h=1}^{g} [1-h/(g+1)]\Big(\sum_{t=h+1}^{n} \hat a_t \hat a{t-h}\Big) \]

e \(ep(\hat \beta_1)=[SE/\hat \sigma]^2\sqrt(\hat v)\), em que SE é o erro padrão usual do MQO para \(\hat \beta_j\).

Desta forma obtêm-se os erros-padrões robustos dos parâmetros. Fazendo uso do pacote sandwich, e a opção para a matriz de variância-covariância consistente com heterocedasticidade e autocorrelação vcovHAC conforme o script abaixo. Conforme a ordem da autocorelação (se acima da primeira ordem), pode ser mais indicado considerar o tratamento para a estrutura geral de Newey e West (1987) e Wooldridge (1989), citados por Wooldridge (2016, seção 12.5).


#Correcao do problema da Autocorrelacao - Duas formas

regressao1 <- lm(log(c) ~log(y) + inf, data=dados)
coeftest(regressao1, vcov=NeweyWest)
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.26490    0.17532    1.51    0.140    
## log(y)       0.95155    0.01844   51.60   <2e-16 ***
## inf         -0.00187    0.00101   -1.86    0.072 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(regressao1, vcov=vcovHAC)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.264900   0.138872    1.91    0.065 .  
## log(y)       0.951550   0.014772   64.42   <2e-16 ***
## inf         -0.001872   0.000928   -2.02    0.052 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


4 Viés de Especificação


Omitir uma variável relevante - Conseqüências:

  1. se a variável omitida for correlacionada com alguma incluída, os betas estimados são tendenciosos e inconsistentes;

  2. se não existir correlação, o intercepto é viesado;

  3. a variância do erro estará errada;

  4. IC, erros padrão e testes de hipóteses não são calculados corretamente.

Inserir uma variável irrelevante - Conseqüências: a perda da eficiência, ou seja, erros-padrão maiores do que as dos parâmetros do modelo verdadeiro.

A conclusão é que incluir uma variável irrelevante é menos problemático do que omitir uma variável relevante. Existem alguns testes para identificar problemas na especificação dos modelos.


4.1 Viés de Especificação: como detectar?


O Teste Reset de Ramsey é usado para analisar má especificação do modelo. É um teste mais geral, tanto de variáveis irrelevantes, quanto de forma funcional incorreta e correlação entre os regressores e o termo de erro. A hipótese nula é de que o modelo inicialmente estimado está bem especificado. Uma falha no teste é que diz que o modelo não está bem especificado, mas não se sabe qual alternativa é a melhor. Deve-se testar outras formas, com outras variáveis e ir refazendo o teste.

Considere um modelo de regressão

\[ y=\beta_0+\beta_1x_1+ \dots +\beta_kx_k+u \]

se a hipótese \(E(u|x_1,x_2, \dots ,x_k)=0\) for violada isto indica que a relação funcional entre as variáveis explicadas e explicativas está mal especificada, como mostra Wooldridge (2018, p. 90).

O Teste Reset adiciona polinômios aos valores estimados por MQO na equação acima para verificar se são significativos e detectar tipos gerais de má especificação de formas funcionais.

Para fazer o teste é preciso definir quantas funções dos valores estimados serão incluídos na regressão expandida. No geral, quadráticos e cúbicos tem dado bons resultados.

A estatistica do teste é dada por

\[ F= \frac{\frac{R^2_{novo}-R^2_{velho}}{m}}{\frac{1-R^2_{novo}}{n-p}} \]

em que m é o número de novos regressores e p é o número de parâmetros no novo modelo.

\(F \sim F_{m,n-p}\) com n igual ao número de observações.


#Regressao Multipla

regressao1 <- lm(lwage ~ educ + I(educ^2) + exper + I(exper^2) + tenure)
summary(regressao1)
## 
## Call:
## lm(formula = lwage ~ educ + I(educ^2) + exper + I(exper^2) + 
##     tenure)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.848 -0.232  0.014  0.251  1.292 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.066469   0.671795    6.05  2.4e-09 ***
## educ         0.289380   0.095023    3.05   0.0024 ** 
## I(educ^2)   -0.007467   0.003321   -2.25   0.0249 *  
## exper        0.005706   0.015446    0.37   0.7120    
## I(exper^2)   0.000596   0.000669    0.89   0.3733    
## tenure       0.009768   0.003091    3.16   0.0016 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.378 on 657 degrees of freedom
## Multiple R-squared:  0.166,  Adjusted R-squared:  0.159 
## F-statistic: 26.1 on 5 and 657 DF,  p-value: <2e-16
#Teste de especificação do Modelos
#Hipotese Nula: Modelo esta bem especificado.
resettest(regressao1, power=2 , type='fitted')
## 
##  RESET test
## 
## data:  regressao1
## RESET = 0.33, df1 = 1, df2 = 656, p-value = 0.6
resettest(regressao1, power=3 , type='fitted')
## 
##  RESET test
## 
## data:  regressao1
## RESET = 1.9, df1 = 1, df2 = 656, p-value = 0.2
#Estimacao do Modelo Linear
regressao2 <- lm(lwage ~ educ + exper + tenure)
summary(regressao2)
## 
## Call:
## lm(formula = lwage ~ educ + exper + tenure)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8574 -0.2360  0.0171  0.2460  1.3081 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.50223    0.12825   42.90  < 2e-16 ***
## educ         0.07572    0.00745   10.17  < 2e-16 ***
## exper        0.01799    0.00408    4.41  1.2e-05 ***
## tenure       0.00986    0.00307    3.22   0.0014 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.379 on 659 degrees of freedom
## Multiple R-squared:  0.159,  Adjusted R-squared:  0.155 
## F-statistic: 41.4 on 3 and 659 DF,  p-value: <2e-16
#Teste de especificacao do Modelos
#Hipotese Nula: Modelo esta bem especificado.
resettest(regressao2, power=2 , type='fitted')
## 
##  RESET test
## 
## data:  regressao2
## RESET = 0.1, df1 = 1, df2 = 658, p-value = 0.7
resettest(regressao2, power=3 , type='fitted')
## 
##  RESET test
## 
## data:  regressao2
## RESET = 0.95, df1 = 1, df2 = 658, p-value = 0.3


O teste Reset de Ramsey é usado para analisar má especificação do modelo. É um teste mais geral, tanto de variáveis irrelevantes, quanto de forma funcional incorreta e correlação entre os regressores e o termo de erro.